﻿using System;

namespace MathService.DSP
{
    [Copyright("Расчет аналогового нормированного ФНЧ Чебышева первого рода", url = "http://www.dsplib.ru/content/filters/ch4/ch4.html")]
    public class ChebyshevLowpassFilterI : CBIH
    {
        public static void Test()
        {
            var wp = 1d;
            var ws = 2d;
            var Rp = 1.5;
            var Rs = 35d;
            var dt = 2d;

            var filter = new ChebyshevLowpassFilterI(wp, ws, Rp, Rs, dt);
        }

        private static double arsh(double x) { return Math.Log(x + Math.Sqrt(x * x + 1)); }
        private static double arch(double x) { return Math.Log(x + Math.Sqrt(x * x - 1)); }

        private static PQ Calculate(double wp, double ws, double Rp, double Rs)
        {
            Func<double, double> eps = R => Math.Sqrt(Math.Pow(10, R / 10) - 1);
            var eps_p = eps(Rp);
            var eps_s = eps(Rs);

            Func<double, double> G = e => 1 / Math.Sqrt(1 + e * e);
            //var Gs = G(eps_s);
            var Gp = G(eps_p);

            var Nd = arch(eps_s / eps_p) / arch(ws / wp);
            var N = (int)Math.Round(Nd);
            if(N < Nd) N++;

            const int r = 0;
            var L = (N - r) / 2;
            var b = arsh(1 / eps_p) / N;
            Func<int, double> A = k => (2 * k - 1) / (2 * (double)N) * Math.PI;
            var a = new double[L].Initialize(i => A(i + 1));

            Func<int, double> Sigma = k => -Math.Sin(a[k - 1]) * Math.Sinh(b);
            Func<int, double> Omega = k => Math.Cos(a[k - 1]) * Math.Cosh(b);

            var sigma = new double[L].Initialize(i => Sigma(i + 1));
            var omega = new double[L].Initialize(i => Omega(i + 1));

            var p = Gp;
            var HQ = new Polynom(1);
            for(var i = 0; i < L; i++)
            {
                var wi = omega[i];
                var wi2 = wi * wi;
                var si = sigma[i];
                var si2 = si * si;
                var w2s2 = si2 + wi2;
                p *= w2s2;
                HQ *= new Polynom(w2s2, -2 * si, 1);
            }
            var HP = new Polynom(p);

            return new PQ(HP, HQ);
        }

        public ChebyshevLowpassFilterI(double wp, double ws, double Rp, double Rs, double dt)
            : this(Calculate(wp, ws, Rp, Rs).P2Z_BilinearTransform(dt))
        {

        }

        private ChebyshevLowpassFilterI(PQ H) : base(H.Q.Coefficients, H.P.Coefficients) { }
    }

    [Copyright("Расчет аналогового нормированного ФНЧ Чебышева второго рода", url = "http://www.dsplib.ru/content/filters/ch5/ch5.html")]
    public class ChebyshevLowpassFilterII : CBIH
    {
        private static PQ Calculate(double wp, double ws, double Rp, double Rs)
        {
            return null;
        }

        public ChebyshevLowpassFilterII(double wp, double ws, double Rp, double Rs)
            : this(Calculate(wp, ws, Rp, Rs))
        {

        }

        private ChebyshevLowpassFilterII(PQ H) : base(H.Q.Coefficients, H.P.Coefficients) { }
    }
}